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Abstract. A novel integration method for quadratic vector fields was introduced by Kahan 
in 1993. Subsequently, it was shown that Kahan’s method preserves a (modified) measure 
and energy when applied to quadratic Hamiltonian vector fields. Here we generalize Kahan’s 
method to cubic resp. higher degree polynomial vector fields and show that the resulting 
discretization also preserves modified versions of the measure and energy when applied to 
cubic resp. higher degree polynomial Hamiltonian vector fields. 


1. Introduction: Kahan’s method for quadratic vector fields 

The study of ordinary differential equations (ODEs) goes back centuries, to the time of Newton, 
Bernoulli, Euler, and contemporaries. Since the invention of the computer in the 1940s much 
attention has been devoted to the best ways to discretize differential equations so that they can 
be solved numerically. Initially, the main emphasis was on all-purpose methods (defined for all 
ODEs), such as Runge-Kutta methods and linear multistep methods, and their quantitative 
accuracy. During the last two or three decades, however, interest has expanded to considering 
special classes of ODEs and purpose-built algorithms that preserve the special features of each 
class. These novel methods are not just quantitatively , but also qualitatively accurate. This 
has resulted in methods that preserve symmetries, first integrals, symplectic structure, measure, 
foliations, Lyapunov functions, etc. These methods are called geometric integration methods 

mm- 

In 1993, Kahan introduced a numerical integration method for quadratic differential 
equations. For the quadratic ODE 

x = f(x) := Q(x, x) + Bx + c, iel", (1) 

(where x £ R ra , Q is an R"-valued symmetric bilinear form, B £ R™ xn , and c £ R") it is defined 
by the map 

x i y x': —-— = Q(x, x') + + x') + c (2) 

where h is the time step. The method © was introduced in [TJ for two examples, a scalar 
Riccati equation and a 2-dimensional Lotka-Volterra system and written down in the general 
form © in [45] (see also references therein). Kahan wrote in the prologue to [14], 



Discretization of polynomial vector fields by polarization 


2 


“I have used these unconventional methods for 24 years without quite understanding 
why they work so well as they do, when they work. That is why I pray that some reader 
of these notes will some day explain the methods’ behavior to me better than I can, and 
perhaps improve them. ” 

Initially, the mystery only deepened, for the Kahan method did not at first sight fit into 
any of the standard methods of discretizing ODEs, nor into any of the new methods that were 
developed as the field of geometric numerical integration grew. Yet in some sense Kahan’s 
prayer has been fulfilled. The Kahan method has been found to have remarkable geometric 
properties. Studies have shown that the Kahan method preserves complete integrability in many 
cases [H 03 EU EH EH Ell EH ED E2J- For a large class of Hamiltonian systems the method 
has a conserved quantity (related to energy) and an invariant measure. It is the restriction of a 
Runge-Kutta method to quadratic vector fields [ 4 ]. However, so far only a part of the observed 
behavior of the method has been accounted for, and the ‘explanations’ to a degree only raise new 
questions, for they reveal aspects of Runge-Kutta methods and of discrete integrability that were 
previously unknown and unsuspected. Maps derived from the Kahan method are birational, with 
birational inverses; thus they are elements of the Cremona group of birational automorphisms. 
The algebra, geometry, and dynamics of this group have been studied extensively 01231, although 
the phenomena illustrated by the Kahan method are apparently new. 

Just one of the unusual features of Kahan’s method is that the formulation (ED is defined 
only for quadratic differential equations. Although its Runge-Kutta formulation is defined for 
all ODEs, the special geometric properties appear to hold only in the quadratic case. Yet there 
is no apparent structure to the set of quadratic differential equations that would distinguish 
them in this way, especially in relation to the birational maps. In this paper we propose 
a natural generalization of Kahan’s method to polynomial vector fields of higher degree and 
show that it does inherit some of the geometric properties— invariant measures, first integrals, 
and integrability—of Kahan’s method in some cases. (An alternative generalization of Kahan’s 
method to higher degree vector fields is considered in m-) 

We first observe that a homogeneous quadratic vector field f(x) can be expressed in terms 
of a bilinear form Q(x,x), as in (ED, using the technique of polarization: 

Q(x, x') = ^ (/( x + x') - f(x) - f{x')) . ( 3 ) 

Then the Kahan method can be obtained by polarizing the quadratic terms of the ODE, 
evaluating them at (x,x'), and by replacing the linear and constant terms by the midpoint 
approximation. 

Polarization is a map from a homogeneous polynomial to a symmetric multilinear form in 
more variables. For example, the polarization of the cubic f{x) is the trilinear form 

1 () i) () 

F( X!,X 2 ,X 3 ) = + ^ X2 + ^ 3 ® 3 )|a =0 

where x, x\, x?, and X3 are all vectors in R". This is equal to A of the coefficient of A1A2A3 in 
f(AiXi + A2X2 + A3X3). It satisfies 

F(x,x,x) = f(x). 

For example, consider the case x G R 3 and write x = (y , z, w) T . Then the polarization of 3 y 2 z is 
2/12/223 + 2/22/321 + 2/32/122 and the polarization of 6 yzw is 2/122W3 + 2/223W1 + V3Ziw 2 + yiz 3 w 2 + 
2/3 Z2W1 + 2/221W3. Polarization was used in [B to obtain linearly implicit, integral-preserving 
methods for Hamiltonian PDEs. 

Polarization of a homogeneous vector field of degree k + 1 will lead to a multilinear form 
in k + 1 variables. We will call these variables Xg, ■ ■ ■ ,Xk, where Xk G R n . The generalization 
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of the Kalian method that we consider in this paper is to evaluate this multilinear form at 
k + 1 consecutive time steps, leading to a fc-step numerical integrator. In this way, the bilinear 
character of the Kahan method carries over to higher degrees. The treatment of the linear 
term x is no longer unique; here we consider the simplest possible option of discretizing x by 
(x k - x 0 )/(kh). 

Definition 1 . Let V = K™ and let F be the multilinear map from V k+1 to K" associated with 
the homogeneous polynomial differential equation 

x = F(x,x,... ,x) (=: f(x)) 

of degree k + 1 on V. The polar map associated with f is the birational map on V k given by 
(xq, ..., x k -i ) (x \,..., x k ) where x k is the solution of the linear system 


Xk Xq 

kh 


F(xq, .. .,x k ). 


( 4 ) 


Note that as both sides of © are linear in Xq and in x k , Eq. ©> the expressions for both x k 
as a function of xq, ..., x k -i, and for xq as a function of x \,..., x k , are rational functions. Thus, 
like the Kahan map, the polar map is birational. However, it is expected that the multilinearity 
of © is more special than mere birationality; when k > 1 there are many birational integrators 
formed from / that are not multilinear. The multistep leapfrog method 


X2 - x 0 
2 h 


= /0 i) 


is an example; maps of this form are not expected to have special geometric properties. 


Proposition 1 . The polar map of a homogeneous quadratic is its Kahan map. If a 
nonhomogeneous quadratic is suspended to a homogeneous form in one dimension higher (e.g. 
if x = x 2 + bx + c is replaced by x = x 2 + bxy + cy 2 , y = 0), then the polarization of the 
suspended vector field, projected to the original phase space, is exactly the Kahan map of the 
nonhomogeneous quadratic. 

Proposition 2 . The polar map is (i) self-adjoint (in the sense of symmetric multistep methods 
m), and (ii) a general linear method restricted to vector fields that are homogeneous polynomials 
of degree k + 1 . 


Proof, (i) Eq. © is invariant under (xo, ■ ■ ■, x k , h) n- ( x k , ■ • •, Xo , —h). 

(ii) From a standard identity in algebraic polarization © p. 110 ], 

F(x 0 ,...,x k ) = ]_ (~l) k+1 ~ m f{x il +... + x im ), 

^ ' ' l<m<fc + l 


where the sum is over all nonempty subsets of { 0 ,..., k}. Using homogeneity of /, we get 


F{x 0: ...,x k ) 


1 

(k + 1 )! 


Y (-1 ) k+1 - m m k+l f 

l<m<fc + l 



m 



There are 2 k+1 — 1 nonempty subsets of { 0 , so in this form, the vector field f(x) is 

evaluated at 2 k+1 — 1 points, each of which is a convex combination of the Xj. These points may 
be taken to be the stage values of a general linear method. □ 
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General linear methods are a natural class of methods that include both Runge-Kutta and 
linear multistep methods [3]- In this case the method has fc — 1 ‘auxiliary variables’ xo, ■ ■ ■, Xk -2 
that are carried forward along with the ‘current point’ Xk-i, but is ‘mono implicit’ in the sense 

that only a single variable, Xk, enters nonlinearly. (When / is degree k + 1, Xk even enters 

linearly.) 

For example, if f(x) is a homogeneous cubic then we can write 

F(x 0 ,x 1 ,x 2 ) = 7 (/(x 0 +xi + x 2 ) - f{x o + X\) - f{x o + x 2 ) - f(x i +x 2 ) 

6 

+ f(xo) + f{x\) + f{x 2 )) 

_ 27 r f Xo + x\ + x 2 \ 8 f f xo + Xi \ 8 I x 0 + x 2 \ 8 f xi+x 2 \ 

- ~& } V 3 ) 6 ; V 2 ) ~ f> ^ 2 ) 6 ; ^ 2 ) 

+ ^f(x o) + ^f{x l) + ~ f (x 2 ). 

The special behavior of the Kahan method is seen most easily on the scalar ODE x = x 2 , 
for which it yields the map xq i —> x\ = Xo/(l — hx o), a Mobius transformation which is easily 
integrated. It can be seen to converge past the singularity at t = l/x(0). In contrast, an explicit 
method (like forward Euler) has no singularity, and an implicit method (like backward Euler) 
does not define a smooth map ip: X -4- X for any sensible domain X C R. We first study the 
polar map associated to a higher-degree analog of this ODE. 

Proposition 3. Let k be a positive integer. The polar map of x = x k+ ^, x £ M, is explicitly 
integrable. 

Proof. Eq. 0 written at time step n becomes in this case 

Xn+k = X n + hk X n X n +i . . . X n +k. 

Dividing both sides by x n ... x n +k, 


X n • • • X n -\-k —1 *^n+1 • • • *£n+fc 

Thus, I n := 1/ (x n ... x n+k - 1 ) obeys 
I n — I n —i hk 

with solution 

I n = Io — nhk. 

Taking logs, log(a: r i) obeys the linear, constant-coefficient, nonautonomous difference equation 
log(x„) + ... + log(x n+fc _i) = - log(/„) 

which is easily solved. □ 

This encouraging behaviour motivates our study of the polar map. Although not exhaustive, 
one large class of vector fields for which the Kahan map is known to have special properties is that 
of the Hamiltonian vector fields. For k = 1 the polar (Kahan) map derived from Hamiltonian 
vector fields on Poisson spaces with constant Poisson structure is known to have a conserved 
quantity and an invariant measure [4] . This explains their integrability in some (low-dimensional) 
cases. We will now show how these phenomena generalize to the case k > 1. 
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We will now consider the case when f(x) is a homogeneous Hamiltonian vector field defined on 
a symplectic vector space. The following result establishes the existence of several first integrals 
of an iterate of the associated polar map. These integrals correspond to ‘modified energies’ of 
the map as they all approximate the Hamiltonian in the limit of small step size. 

Proposition 4. Let H.V —>• R, the Hamiltonian, be a homogeneous polynomial of degree k + 2 
on R” and let H be a symmetric (k + 2)-tensor such that H(x) = -jj^yH^XjX,... ,x). Let H be 
a constant invertible antisymmetric n x n matrix , and let uj be its associated symplectic form on 
V, i.e., 

uj: V x V —> R, uj(u , v) = u T flv. 

Let K = fl _1 . Then 

(i) the Hamiltonian ODE on V associated with (H,fl) is 

±= + ( 5 ) 

(ii) the associated polar map is defined via 

xjL ur = why KHix ° . (6) 

(Hi) the associated polar map has k independent k-integrals 

oj(xq, X\), uj(x 1 ,x 2 ),..., uj(x k -i,x k ). (7) 

Proof. First note that the components of the gradient VH(x) are homogeneous polynomials of 
degree k + 1 in the components of x and we have 

(jfcT2j! H^’ *’-’*) = = kh ViJ(x)T ’ 

SO 

(x,...,x,-) = VH(x). 

This yields (i) and (ii). 

Recall that a ^-integral of a map is an invariant of the fcth iterate of the map [5]. We first 
show that uj(xo,x\) = uj(x k ,Xk+i)- Writing a = kh/(k + 1)!, we have 

lo(x 0 ,Xi) - uj(x k ,x k +i) = ui(x k - aKH(x 0 ,... ,Xk,-),x k +i ~ aK H(xi,..., x k +i, •)) ~ ui(x k ,Xk+i) 
= -auj(K H(x 0 ,.. .,Xk, -),Xk+ 1 ) + aw{KH(x\,.. .,x k +i, -),x k ) + 
a 2 ui(K H(a.’ 0 ,.. .,x k , -),K H(xi,... ,x k +i, •)) 

— aH (^0; * - * i Xk, Xk -\-1 ) ttH (x± , . . . , X k -{- 1, Xk ) 

+ a 2 H(a;i,..., x k +i, KH(x 0 , ■■■,x k , •)) 

— aH (xfl ,..., Xk, Xk~k i ) aH (x \,..., Xk+i > k) 

+ a 2 H(zi .. .,x k+ i,(x k - x 0 )/a) 

= 0 . 

Now if I{x) is any /c-integral of a map <p, then so is I o (p( m > for any integer m [5]. This yields 
the remaining fc-integrals and concludes the proof. □ 
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The integrals all approach the Hamiltonian of the original system as x m —> x(mh ), for in 
this limit we have 


1 


h(k + 2) 


^ Om ? %m+ 1) — 


1 


7^{Xrm%m +1 X rn ) 


—>■ 


h(k + 2) 

^ ^ uj(x(mh),x(rnh )) 

V-ff (x(mh)) T x(mh) 


k + 2 
= H(x(mh)). 


Note that the /c-integral uj(xk-i, Xk) is a rational function of 0o> • • ■ , Xk-i), for Xk is defined 
through ©; the other fc-integrals are all quadratic. 

Analogous results hold for Poisson systems of the form ([5]) where K is antisymmetric but 
not invertible. These can be established by performing a linear change of variables that puts K 
in its Darboux normal form. Because the discretization method is linear, it commutes with this 
change of variables. 


3. Invariant measure of the polar map 


The next proposition shows that the polar map in this case is also measure-preserving in the 
extended phase space. The measure is ‘modified’ in the sense that it converges to the invariant 
measure of the ODE in the limit of small step size. 

Proposition 5. Let K be a constant antisymmetric n x n matrix and let H:P fe+2 —> R be 
multilinear. Let p be a constant measure on V and let p k be the corresponding product measure 
on V k . Then the map on V k induced by the polar map 0 associated with the homogeneous 
Hamiltonian vector field x = K H(x, a;,..., x, ■) has the invariant measure 

,k 




where 


det(I - cK H (xo, ■ ■ -,x k -i , •)) 
h 


c = 


(*-!)!' 
Proof. First note that we have 
1 


(fc + 2)! 


H(a 


) = H(x) = 


1 


k + 2 


VH(x) t x = 


1 


(k + 2 ) 0 + 1 ) 


n H"{x) 3 


and 


So 


VH(x) = - - H"(x)x. 

w jfc + 1 w 

1 


0 + 1 )! 


H(+..., + •) = VH( x), 


k\ 


H(x,...,x,-,-) = H"{x). 


Let X = [xf ...., x'^_ 1 ] T . We want to prove that the Jacobian of the map 

p: Oo,..., x k - 1 ) (x' 0 ,..., x'k-x) 


= x i 
= X 2 


( 8 ) 


f4_! = Xk = X 0 + p^A'HOo,... ,x k , •) 
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det 


dip _ det(I - cK H(si,. ..,x k ,•)) 
d X det (I - cKH(x 0 , ■■■, x k -i, •, •)) 


We first observe that 


dip , dx' 

del M =de ‘^„ 




dx n 


This follows directly from the format of the Jacobian and in fact 


i dp 

det — = det 
aX 


O 

O 


o 

dx k 

dxo 


I 

O 


o 

dock 

dx\ 


o 

I 


o 

dock 

dx k -2 


o 

o 


I 

dx k 

dx k - 1 


= det 


dx k 

dx n 


Differentiating (0 on both sides with respect to xq 1 and using the symmetry of H we have 


dx k 
dx o 


dx]$ 

I + cKH{x 0 ,... ,x k -i,-,-) ~z. -h cKH(xi,... ,x k , •, •) 

OX 0 


dx o 
dx 0 


Rearranging the terms we obtain 


dx k 

dx 0 


(/ - cK H(x 0 ,..., a:fc_i, •)) 1 {I + cKH(x 1 ,..., x k , ■)), 


and 


det - det (-f + cK H(xi, ...,x k ,-, •)) 
dx o det(/ - c/v H(a;o,... ,x k - i, ■, •))' 


Using det(H) = det(H T ) and the Sylvester determinant theorem det(J + AB) = det(J + BA) in 
the numerator, we obtain 


det - det (^ ~ cK H(xi, ...,x k ,-, •)) 
dx o det(/ - cA' H(x 0 ,... ,x k -i,-, •)) ’ 

establishing the result. □ 

Note that in the case k = 1, in which case 0 reduces to the Kahan method for homogeneous 
cubic Hamiltonians, the invariant measure 0 can be written as p/ det (I — A hf(x )), which is 
the form of the invariant measure for the Kahan method found in j?j. 


4. Integrability of the polar map 

The next property of the polar map concerns a {k — l)-dimensional symmetry group, so it is a 
phenomenon that only appears for k > 1. 

Proposition 6. (i) The kth iterate of the polar map 0 is equivariant with respect to the 

scaling symmetry group x m X m x m , m = 0,..., k — 1, where A m = 1, i.e., the map 

0 has a (k — 1)-dimensional k-symmetry group. 

(ii) The measure (B is invariant under this scaling group. 

(Hi) The integral 11^=0 u ( x m, Xm+i) is invariant under this scaling group. 
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(iv) When k is even, the 2-integrals 

ljj{x 0 ,X 1 )u}(X 2 ,X 3 ) . ..U}(Xk- 2 ,Xk-l) 


and 


U}(x 1 ,x 2 )uj(x3,x i ) ■ ■ ■ u(Xk-l,X k ) 

are invariant under this scaling group. 


Proof, (i) Under the map x m K > A mX m , the final equation defining the polar map, ©■ is 
transformed to 


Aq Xk = Aq^o + khK 


(k + l)\ 


H(Aq^0: Ai3"i, ■ • ■ , Ak—lXk—1, A QXk, *) 


which is identical to © under the condition J| m _ 0 A m = 1. Therefore, the function 
tp(xo, ..., Xk-i) (= Xk) defined through the solution of © scales as <p(xo, ■ ■ ■, Xk-i) 
Ao<p(xo, ■ ■ ■, Xk- 1 ). The fcth iterate of the polar map can be written 


xo k) = (p(x 0 ,. ■ .,x k -i) 
x[ k) = (p(x 1 , . ..,Xk) 


a4_l = <p(Xk- 1, • • -,X2 k-2) 

and each equation is invariant under the action x m >->• A mo d(m,k) x m induced on the iterates 
of the map. 

(ii) Follows from H(xo,.. .,x k -i) = H(A 0 a;o,..., Xk-ix k -i). 

(iii) We have 

k— l k— l 

II UJ ( X ™,Xm+ 1) ^ n cj ( A m x m , A m+ i£ m+ i) 

m= 0 m— 0 

k -1 k -1 

— j| A m 11 *^m+l) 

m=0 m—0 
fc-1 

= w(x m ,x m+ i) 
m—0 

which establishes the result. 

(iv) Under the symmetry, each of the given 2-integrals is multiplied by a factor rim^o w hi c h 
establishes the result. 

□ 


These results yield a 5-parameter family of integrable 4-dimensional rational maps. 
Corollary 7. The polar map is completely integrable in the case k = 2, n = 2. 

Proof. The 2nd iterate of the polar map in this case has a 1-dimensional measure-preserving 
symmetry group. The map thus descends to a measure-preserving map on the 3-dimensional 
quotient na. The two integrals of the 2nd iterate of the polar map are invariant under the 
symmetry and hence also pass to the quotient. This yields a 3-dimensional measure-preserving 
map with 2 integrals, thus integrable [2]. The reconstruction dynamics obey a 1-dimensional, 
linear, constant-coefficient, nonautonomous difference equation and hence are integrable. From 
the integration of the 2nd iterate, the integration of the polar map itself is immediate. □ 
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Example 1. To be fully explicit we give here the integrable rational map obtained in the case 
k = 2, n = 2. Let ( q,p ) be coordinates on V = R 2 and let the Poisson tensor and Hamiltonian 
be 

' 0 r 

-1 0 

Then the polar map on V 2 is (qo,Po,qi,Pi) | -t (<7i, Pi, < 72 , P 2 ), with 

<7o + 2h (bq^q ± + c(2p 0 <7o<7i +Piq%) + d(2p 0 p 1 q 0 + plq x ) + eplpi) 


K = 


H = aq 4 + Abq 3 p + 6 cq 2 p 2 + Adqp 3 + ep. 


<72 = 

P2 = 


1 - Ah 2 A 

Po - 2 h (aqlqi + b(2p 0 q 0 q 1 + Piqp + c(2p 0 pig 0 + pggi) + dppPi) 

1 - Ah 2 A 


where 


A = 

c d 
d e 

2 2 , 
P 0 P 1 + 

b c 
d e 

(PoPiqi +p 0 plqo) + 

b c 
c d 



+ 

a b 
c d 

(Pi<7o<7i +Po<7o<7i) + 

a c 

c e 


(Po<h + 7h<7o) 


PoPiqoQi 

The map is birational of degree 3 over degree 4. The two 2-integrals are 


(9) 


2 2 
<7o<7i- 


qoPi - qiPo and q ± p 2 -P 192 , 

where q 2 and p 2 are given in ©. The invariant measure is 
dqo A dpo A dq± A dpi 
1 - Ah 2 A ‘ 

If the degree k > 2 or the dimension n > 2 then the geometric properties described above 
are not enough to ensure integrability. Indeed, we find that the polar map associated with a 
homogeneous planar quintic Hamiltonian (i.e. k = 3, n = 2) does not pass the entropy test for 
complete integrability mm- 
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